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We extend our previous results of the 14th post-Newtonian (PN) order expansion of 
gravitational waves for a test particle in circular orbits around a Schwarzschild black hole 
to the 22PN order, i.e. v 4,4 beyond the leading Newtonian approximation where v is the 
orbital velocity of a test particle. Comparing our 22PN formula for the energy flux with 
high precision numerical results, we find that the relative error of the 22PN flux at the 
innermost stable circular orbit is about 10 -5 . We also estimate the phase difference between 
the 22PN waveforms and numerical waveforms after a two-year inspiral. We find that the 
dephase is about 1CT 9 for fj,/M = 1CT 4 and 1(T 2 for fi/M = 1CT 5 where fj, is the mass of the 
O ■ compact object and M the mass of the central supermassive black hole. Finally, we construct 

a hybrid formula of the energy flux by supplementing the 4PN formula of the energy flux 
for circular and equatorial orbits around a Kerr black hole with all the present 22PN terms 
^jT^ for the case of a Schwarzschild black hole. Comparing the hybrid formula with the the full 

numerical results, we examine the performance of the hybrid formula for the case of Kerr 
black hole. 

^ . 81. Introduction 

m 

m 

■ Gravitational waves (GWs) are expected to be detected in this decade by up- 

*n . coming ground based detectors such as Advanced LIGO, 1 ) VIRGO 2 ) and KAGRA. 3 ) 

One of the most promising sources of the GWs for the ground based detectors is 
the coalescing stellar-mass compact object binary system. For the detection of the 
GWs and the subsequent parameter estimation of the source, one will correlate the 
noisy signal of the detector with the template bank of theoretical waveforms. Thus, 
it is important to predict the waveforms very accurately and efficiently not only for 
developing the data analysis strategy but also for extracting the physical information 
\ of the source. 

A conventional method to predict inspiral waveforms from coalescing binaries is 
the post-Newtonian (PN) expansion of the Einstein equations, 4 ) in which the relative 
velocity of the binary v is assumed to be smaller than the speed of light. Currently, 
the amplitude (orbital phase) of gravitational waves is derived up to 3PN (3.5PN), 
i.e. v 6 (v 7 ) beyond the leading order, 5 )~ n ) for the non-spinning compact binaries in 
quasi-circular orbits. (Note that the 3.5PN amplitude for I = m = 2 mode is derived 
in Ref. 12).) Although the PN expansion accurately describes the early phase of the 
inspiral, the approximation breaks down in the late phase of the inspiral. For the late 
inspiral and the subsequent merger and ringdown phases one needs to calculate the 
waveforms using a full numerical solution of the Einstein equations. 13 )' Since the 
computational cost to perform full numerical simulations which include the whole 
process of the coalescence is too expensive, one should make theoretical templates by 
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matching the PN waveforms in the early inspiral phase with the numerical waveforms 
in the late inspiral and the subsequent merger and ringdown phases. 15 ^' 16 -* Thus, from 
the point of view of the computational cost it is important to obtain higher PN order 
expressions of the GWs and to investigate whether the high PN order expressions 
extend the region where the PN expansions of GWs are sufficiently reliable. 

For this purpose, we derive the high PN order expressions of the GWs focusing on 
a binary system which consists of a compact object of mass (i in circular orbits around 
a Schwarzschild black hole of mass M by assuming fi <C M. Such extreme mass ratio 
inspirals (EMRIs) are one of the main candidates of gravitational waves for space- 
based detector such as eLISA. 17 ) Another approach to model EMRIs is the black 
hole perturbation formalism, which uses the mass ratio as an expansion parameter. 
Although the black hole perturbation theory is limited to the case \x <C M, it is rather 
easier to compute the high PN order expansions of the GWs than using the standard 
PN approximation. Moreover, one can use the numerical results of the black hole 
perturbation theory to investigate the relative accuracy of the PN expansions since 
in the black hole perturbation theory there are no assumptions on the orbital velocity 
of the compact object. 

Using the first order black hole perturbation theory, the gravitational waveforms 
and energy flux to infinity for a test particle in circular orbits around a Schwarzschild 
black hole were computed up to 1.5PN in Ref. 18). (See Refs. 19)-21) for recent 
reviews for orbital evolution due to the small body's interaction with its own grav- 
itational field, i.e. gravitational self-force.) The calculation of the energy flux is 
extended up to 2.5PN by numerical fitting in Ref. 22). Then, in Ref. 23) the 4PN 
expressions of the energy flux were derived by fitting with a very accurate numer- 
ical calculation of the energy flux. It was also found that logf terms appear at 
3PN and 4PN in Ref. 23). These logv terms were confirmed in Ref. 24), which 
computed analytically the 4PN expressions of the energy flux. These calculations 
were extended to 5.5PN for the energy flux in Ref. 25) and for the waveforms in 
Ref. 26). Recently the 14PN gravitational waveforms for a test particle in circu- 
lar orbits around a Schwarzschild black hole have been derived. 27 ) By comparing 
the 14PN energy flux with numerical results, it is shown that the relative error of 
the post-Newtonian energy flux becomes smaller when the PN order is higher. It 
is also shown that the 14PN waveforms using factorized resummation suggested in 
Ref. 28) will provide the data analysis performances of EMRIs comparable to the 
ones resulting from high precision numerical waveforms. In this paper we extend the 
14PN results in Ref. 27) to the 22PN, which is the highest order we can derive with 
reasonable time using our current code.*) We find that if one does not use any re- 
summation technique the 22PN gravitational waveforms will be necessary to achieve 
the data analysis accuracies comparable to the ones using high precision numerical 
waveforms. 

*^ We note that the number of terms necessary to derive the PN expressions grows exponentially 
when the PN order becomes higher. Our current code uses 70, 3.3 x 10 2 , 1.9 x 10 3 , 1.1 x 10 4 and 
3.3 x 10 4 MBytes memory, taking seconds to a few days, to compute multipolar waveforms for 
I = m = 2 mode at 6PN, 10PN, 14PN, 18PN and 22PN respectively. Thus, it is difficult to obtain 
23PN or higher order expressions with reasonable time by using our current code. 
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This paper is organized as follows. In Sec. 2, we start with a brief summary of 
the first order black hole perturbation formalism. Sec. 3 describes a brief summary of 
an analytic formalism to compute homogeneous solutions of the Teukolsky equation. 
Sec. 4 is devoted to the comparison of our 22PN expressions with high precision 
numerical results: The total energy flux to infinity is compared in Sec. 4.1, the 
phase of gravitational waveforms during two-year inspiral Sec. 4.2 and the total 
energy flux to infinity for the case of a Kerr black hole Sec. 4.3. We conclude with 
a brief summary in Sec. 5. Since the 22PN expressions are too large to be shown in 
the paper, we only show the 7PN expression of the total energy flux to infinity in 
Appendix A. In Appendix B we show the 22PN expression of the total energy flux 
to infinity, which can be used for numerical computation in double precision. The 
22PN expressions for the modal energy flux to infinity and factorized waveforms 28 ) 
will be publicly available online. 29 ^ Throughout this paper, we work in the units of 
c = G = 1. 

§2. Teukolsky formalism 

We consider the gravitational waves emitted by a test particle moving on circular 
orbits around a Schwarzschild black hole using the Teukolsky formalism. 30 ) In this 
section, we recapitulate necessary equations following the notation in Ref. 31). 

In the Teukolsky formalism, the gravitational perturbation of a Schwarzschild 
black hole is represented by the Weyl scalar ^4, which represents the gravitational 
waves going out to infinity as 

$4 ->■ 7^(h+ - ihx), for r — > 00, (2-1) 

where dot ' represents a time derivative, d/dt. 

The Weyl scalar can be separated into radial and angular parts if we decompose 
<p4 in Fourier harmonic components as 



1 f°° 

-Y, due-^Re^ir) _ 2 Y em (6,<p), (2-2) 

^ n „„ J —OO 



where £ > 2, — i < m < £ and s Yg m (9,(p) is the spin-weighted spherical harmonics. 9 ) 
The radial function Re muJ (r) satisfies the inhomogeneous Teukolsky equation, 



2 d ( 1 d 



dr V Adr 



Rlmu{r) = T tmul {r), (2-3) 



with A = r(r — 2M) and 



r 2 



U(r) = — [u} 2 r 2 - 4iu(r - 3M)] -(£-!)(£ + 2), (2-4) 



where Te mu) is the source term which is constructed from the energy momentum 
tensor of the small particle. 
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We use the Green function method to solve the inhomogeneous Teukolsky equa- 
tion Eq. (2-3). The outgoing- wave solution of the Teukolsky equation Eq. (2-3) at 
infinity is given by 

Rtmco{r -> oo) = - / dr- 



^r 3 e^*Z £mw , (2-5) 

where r* = r + 2Mln(r/2M — 1) and -R^ w is the homogeneous solution of Eq. (2-3) 
which satisfies the ingoing-wave boundary condition at horizon as 

( Dtrans A 2 -iur* f * . _ 

_ I ^imuj ^ e IOr T OO, , . 

^ \ r 3 Bg iu) e iur * + r- l B'™ w e- iwr * for r* -»• +00. 1 J 

When the test particle moves on a circular orbit around a Schwarzschild black 
hole, the angular frequency J7, the specific energy E and the angular momentum L 
of the particle are given by 



M r -2M . = VMF B 

' 4 v"'o( r o - 3M) v/1 - 3A//r 



where ro is the orbital radius. The frequency spectrum of T^ muJ has peaks at the 
harmonics of the orbital frequency 00 = mil. Then one finds that Zi muJ in Eq. (2-5) 
takes the form 

Zirrwj = Z £muJ 5(uj - mtt), (2-8) 

where 



7 I™ 

J (.mu 



-obe m — 2i_ibe m I 1 + 



2 (ro - 2M ) 



euro ( M 1 . 
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+ y -lhm ~-2 hm ( 1 + rQ _ 2M ) f t M 
+ 2-2 i '«m'od<J ( r o) 

and prime, denotes differentiation with respect to r and s bi m are defined by 



ohm =\ [(t ~ + W + 2)] 1/2 oY £m (|, 0) (2-10a) 

= [(£ - 1)(* + 2)] 1 / 2 _!^ m (~ o) -, (2- 10b) 

V2 / r 

=-2^m (J, °) ( 2 - 10c ) 
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Using the amplitudes Zi mw , the gravitational wave luminosity is given by 



alE 
~dt 
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(2-11) 



1=1 m=-i 



Choosing {6,(p) as the angles defining the location of the observer relative to the 
source, the gravitational waveforms can be expressed as 



where uj = mfi is the frequency of the gravitational waves. 

We derive the post-Newtonian expansions of the gravitational wave luminosity 
Eq. (2-11) and gravitational waveforms Eq. (2-12) by computing the amplitude Zi rnw . 
For the computation of Z^ mu , one needs to obtain the series expansion of the ingoing- 
wave Teukolsky function R^ lu in terms of e = 2Mw = 2MmQ = 0(i> 3 ) and z = 
ujr = 0(v) and the asymptotic amplitudes in terms of e, where v = (M/ro) 1 / 2 . 

We use a formalism developed by Mano, Suzuki and Takasugi 32 )' 33 ^ to compute 
Rfrnuj an d Bfy^. In the next section, we give a brief review of the formalism for the 
convenience of the reader. 

§3. Analytic solutions of the homogeneous Teukolsky equation 

We adopt the formalism developed by Mano, Suzuki and Takasugi (MST) 32 )> 33 ) 
to obtain high post-Newtonian order expansion of gravitational waves. In the MST 
formalism, the homogeneous solutions of the Teukolsky equation are expressed using 
hyper geometric functions around the horizon and Coulomb wave functions around 
infinity. Since the matching of the two kinds of solutions can be done analytically 
in the overlapping region of convergence, one can obtain analytic expressions of the 
homogeneous solutions without numerical integration. Moreover, the series expan- 
sions are naturally related to the low frequency expansion (See Sec. 3.2). Thus, the 
formalism is very powerful to compute high post-Newtonian order expansion of grav- 
itational waves. Using the MST formalism, the energy absorption of gravitational 
waves into the horizon induced by a particle in circular orbits around the equatorial 
plane of a Kerr black hole was calculated up to 6.5PN order beyond the quadrupole 
formula. 34 ) The energy flux of gravitational waves to infinity by a particle in slightly 
eccentric and inclined orbits around a Kerr black hole was also computed up to 
2.5PN order in Refs. 35), 36). In Rcfs. 26) and 37), the author of this paper was part 
of collaborations that applied the formalism to obtain gravitational waveforms for a 
test particle in circular orbits around a Schwarzschild black hole up to 5.5PN and a 
Kerr black hole up to 4PN. Extending the formalism to very high post-Newtonian 
order, we derived the 14PN expressions of gravitational waves for a test particle in 
circular orbits around a Schwarzschild black hole in Ref. 27). Although the MST 
formalism can be applied to the case of a Kerr black hole, we set q = 0, where q 
is the non-dimensional spin parameter of the Kerr black hole, since we are dealing 



h + — ih x = 




y £m (0,V9)e^*-<) 



(2-12) 
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with the the case of a Schwarzschild black hole in this paper. We refer the interested 
reader to a review Ref. 31) for details of the formalism. 

3.1. Ingoing-wave solution 

A homogeneous solution of the Teukolsky equation in a series of Coulomb wave 
functions Rq is defined as 

2— °° ( ~\ ' \ 

* = * " l) £ H)"^3^f <CW2i - e,,). (34) 

n=— oo x ' 

where z = wr, (a) n = F(a + n)/r(a), and Fn(j],z) is a Coulomb wave function 
defined by 

Fjvfa, z) = e- iz 2 N z N+1 r{ ^N + 2) ) ^ N + l ~^ 2N + 2 5 2 «), (3-2) 

where ^(a, /?; z) is a confluent hyper geometric function. 38 ) Observe that the so-called 
renormalized angular momentum v is introduced in the homogeneous solution in a 
series of Coulomb wave functions Eq. (3-1). This parameter v is a generalization 
of I and v — > I as e — > 0. The parameter is determined by the conditions that the 
series converges and actually represents a homogeneous solution of the Teukolsky 
equation. 

Substituting the solution in a series of Coulomb wave functions Eq. (3T) into 
the Teukolsky equation Eq. (2-3) with Ti muJ = 0, we obtain a three-term recurrence 
relation for the expansion coefficients a v n 



«+i + « + 7X-i=0, (3-3) 



where 



v _ ie(n + u — 1 + ie)(n + v — 1 — ie)(n + u + 1 + ie) 
° n ~ (n + i/ + l)(2n + 2i/ + 3) ' ( ' 

F = -£(£ + 1 ) + (n + v)(n + v + l) + 2e 2 + C ( ^ - " \ (3-4b) 

v _ ie{n + v + 2 + ie)(n + v + 2 — ie){n + v — ie) . 
7n " {n + v){2n + 2u- 1) ' ( ° j 

The series Eq. (3T) converges if the renormalized angular momentum v satisfies 
the following equation, 31 )' 32 ) 

RnLn-l = 1, (3-5) 

where R n and L n are defined by continued fractions, which are derived by the three- 
term recurrence relation Eq.(3-3), as 

Rn _ < In _ In "n.7n+l <+l7n+2 ^ 
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(3-7) 



-1 h'n t In^n-L H n 

From the continued fractions, R n and L n , one can obtain two kinds of the expansion 
coefficients, a v n . In general, these two kinds of the expansion coefficients do not 
coincide. Eq. (3-5) is a condition such that the two types of the expansion coefficients 
coincide and the series converges. 31 )' 32 -* If we choose v using Eq. (3-5), the series 
of Coulomb wave function Eq. (3-1) converges for r > r+. From Eq. (3-4), we find 
aC^ -1 = "fn and (3Zn~ l = fin so that gC^ _1 satisfies the same recurrence relation 
Eq. (3-3) as a v n does. This shows that Rq U ~ 1 is also a homogeneous solution of the 
Teukolsky equation, which converges for r > r+. 

As for the homogeneous solution in a series of hyper geometric functions, which 
converges for r < oo, the expansion coefficients satisfy the same three-term recur- 
rence equation (3-3) derived by using a series of Coulomb wave functions. Thus one 
can use the same renormalized angular momentum v to compute both series. This 
fact is important to match the solution in a series of Coulomb wave functions, which 
converges at infinity, with the one in a series of hypergeometric functions, which con- 
verges at the horizon. As a result of the matching, one can obtain the ingoing-wave 



solution Rf^u' which converges in the entire region, as 



R in 



Imu) 



K U R U C + K. 



v-\Rq 



(3-8) 



where 



A" 



= (26)" 



-»- N 2 2 i N r{3 - 2ie)r(N + 2u + 2) 



r(N + v + 3 + ie)r(N + v + 1 + ie)r(N + v-l+ie) 




, r(n + N + 2v + 1) F{n + + ie)r(n + v + 1 + ie) 



(n - 
(-1) 



N)\ 



r(n + u + 3 



ie)r(n + v + I 

x -1 



ie) 



I 



ie, 



(N - n)l(N + 2v + 2) n {v + Z + ie] 



-a" 



(3-9) 



and N can be any integer. The factor K u is a constant to match the solutions in the 
overlap region of convergence and independent of the choice of N . 



One can obtain analytic expressions for the asymptotic amplitudes B 



trans 

Imu) ' ^Imuj 



B) 



and Bj™; defined in Eq. (2-6) by comparing the asymptotic behavior of RfL^. Eq. (2-6) 
with Eq. (3-8) at r* 



ptrans 
Imui 



±oo. They are derived as 



(3-10a) 



where 



Imu) 

rircf 

Imw 



ie 



,sin7r(V + ie) 
sin7r(i/ — ie) 
u 3 [K u + ie^ v K. v ^\AZe itX ^, 



K. 



A': 



-ie ln t 



(3-10b) 
(3- 10c) 
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I> + 3 + ie) 
r(v-l- ie) 



(3-lla) 



8 



Fujita 



ML = 2^e-¥ e -f ^-D £ (-1)" , * " " C (3-llb) 

n=— 00 

We summarize how to compute the gravitational waves using the MST formal- 
ism. One first should determine v by solving Eq. (3-5). Then, one can compute 
the expansion coefficients a y n using the continued fractions Eq. (3-6) for n > and 
Eq. (3-7) for n < with the condition o,q = % u ~ 1 = 1. Substituting the expansion 
coefficients a v n into Eqs (3-8) and (3-10), one can compute the ingoing- wave solution 
of the Teukolsky equation R l g mui and the asymptotic amplitudes Bf^, which are 
necessary ingredients to compute Zg mw in Eq. (2-9) and hence allow one to compute 
the energy flux to infinity Eq. (2-11) and gravitational waveforms Eq. (2-12). 

3.2. Low frequency expansions of solutions 

In the practical calculation to determine v, we solve an alternative equation 
which is equivalent to Eq. (3-5) for n = 1 

^ + 0^1 + 70^-1 = 0, (3-12) 

where R\ and L_i are given by the continued fractions Eq. (3-6) and Eq. (3-7) 
respectively. 

In the followings, we consider the low frequency approximation for Eq. (3-12). 
In the limit of low frequency, we solve Eq. (3-12) by requiring v — > £ for e — > 0. 
Since the orders of a v n) 7^, (3%, R\ and L_i in e are O(e), 0(e), O(l), 0(e) and 
O(e) respectively except for certain values of n < 0, 31 )' 32 ) one can systematically 
compute the low frequency expansion of v. The closed analytic form of v at 0(e 2 ) is 
given in Refs. 31), 32). The expansion coefficients a v n can be computed by using R n 
and £-| n |) whose order in e are O(e) for sufficiently large n. Therefore, the order of 
the expansion coefficients a v n in e are 0(e' n ') for sufficiently large \n\. Using the low 
frequency expansion of a^, one can easily derive that of the asymptotic amplitudes 
Eq. (3-10). Since the homogeneous solution of the Teukolsky equation Rj^^ Eq (3-8) 
is a function of z = 0(f) and ~ 

0( e M) = 0(w 3 l"l), the homo geneous solution is 
naturally related to the post-Newtonian approximation. For the calculation of Rf muj 
in the post-Newtonian approximation, however, one should add a larger number of 
terms in Rq Eq. (3-1) for n < than that for n > 0. This is because the post- 
Newtonian order of a series of Coulomb wave functions for n < 0, (a u _^F_^ +u (2 i — 

e, z))/(aQF u (2i - e,z)) ~ 0(v 2 ^), grows slower than that for n > 0, (a v n F n+u (2i - 
e, z)) / (aQF v (2i — e, z)) ~ 0(v in ) . 26 ) Then, it is straightforward to compute the 
post-Newtonian expansion of the energy flux to infinity Eq. (2-11) and gravitational 
waveforms Eq. (2T2) using the post-Newtonian expansion of Z£ mLJ Eq. (2-9). The 
closed analytic form of Z^ muJ at 2.5PN for the case of a test particle moving in circular 
orbits around a Schwarzschild black hole is given in Ref. 26). 
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§4. Comparison with numerical results 

4.1. Total energy flux to infinity 

Fig. 1 shows the comparison of the energy flux to infinity between a high pre- 
cision numerical calculation and post-Newtonian approximations up to 22PN. The 
numerical computation is based on a technique in Ref. 39), 40). The accuracy of the 
numerical calculation is mainly limited by truncation of ^-mode. In this work we set 
I = 25 which gives relative error better than 10~ 14 up to the innermost stable circu- 
lar orbit (ISCO). n-PN flux needs I up to n + 2. We note that the agreement between 
the numerical energy flux and post-Newtonian energy flux becomes better when the 
PN order is higher even around ISCO. The relative error of the 22PN energy flux 
around ISCO is about 10 -5 and an order of magnitude better than that of 14PN 
energy flux, which was derived in our previous paper. 27 ) We also note that the ac- 
curacy can be improved if we compute the energy flux using factorized resummation 
waveforms. 28 ) If one uses the factorized resummation waveforms, the relative error 
of the energy flux around ISCO can be reduced to 10 -6 for 14PN expressions 27 ) and 
10" 7 for our 22PN expressions. 
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Fig. 1. Absolute values of the difference of the energy flux to infinity between numerical results and 
the PN approximation as a function of the orbital velocity. Note that the relative error of the 
energy flux between 22PN and the numerical calculation is an order of magnitude smaller than 
the one between f 4PN and the numerical calculation even around ISCO, v = = 0.40825. 



4.2. Phase difference during two-year inspiral 

From Fig. 1, one may expect that 22PN expression for gravitational waves can 
be used for data analysis of EMRIs since the relative error with a high precision 
numerical calculation is about 10~ 5 around ISCO. To investigate the applicability of 
the PN expressions to the data analysis, we compute the phase difference during two 
years quasi-circular inspiral between the PN waveforms and the numerical waveforms. 
Following Ref. 41), we choose two systems, denoted as System-I (System-II), which 
correspond to the early (late) inspiral phase of an EMRI in the frequency band of 
eLISA. System-I has masses (M, /i) = (10 5 ,10)Mq and inspirals from ro — 29M 
to ro ~ 16M with associated frequencies /gw G [4 x 10 -3 , 10 _2 ]Hz. System-II has 
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masses (M,fi) = (10 6 , 10) Mq and starts inspiral from tq ~ UM to r$ ~ 6.0M with 
frequencies /gw G [1-8 x 10~ 3 ,4.4 x 10~ 3 ]Hz. System-I (II) has the mass ratio of 
10~ 4 (10 -5 ) and ~ lx 10 6 (~ 5 x 10 5 ) rads of the orbital phase due to the two-year 
inspiral. 

For the calculation of the phase, we use the method described in Ref. 42), which 
is also used in Refs. 26), 41). The phase of the waveforms can be described as 
m Jq fi(t')dt' — <Z^ m (i), where Q{t) = ^/ M/r^it) and Wi m is the phase of Zg m ^. To 
compute the radius of the orbits as a function of time ro(i), we use a interpolation 
method to save computational time. For the interpolation, using total energy flux 
dE/dt induced by a particle we compute dr^/dt for 10 3 points data of v from v = 0.01 
to v = 0.408. Then, from the 10 3 points data of (v, dr^/dt), we compute (ro(t), 
^imit)) using the cubic spline interpolation. 43 ) 

Fig. 2 shows the absolute values of the phase difference of the dominant mode 
/i2,2 between the PN and the numerical calculation during the two-year inspiral. 
For System-I (II), the absolute values of the dephasing between the PN waveforms 
and the numerical waveforms after the two-year inspiral are about 4 x 10 1 (3 x 10 3 ), 
9xl0~ 3 (5X10 1 ), 10^ 5 (1), 8xl0^ 9 (2xl0~ 2 ) and 10" 9 (10" 2 ) rads for 5.5PN, 10PN, 
14PN, 18PN and 22PN respectively. Using the factorized resummation waveforms, 
the absolute values of the phase difference for System-I (II) can be reduced to 3 x 10 -8 
(9x 10 -4 ) rads for 14PN waveforms 27 ) and 10 -9 (10 -4 ) rads for our 22PN waveforms. 

Since System-II represents the inspiral in the most strong-field of a Schwarzschild 
black hole, the dephase between the 22PN waveforms and the numerical waveforms 
after two-year inspirals is expected to be smaller than 10 -2 rads for extreme mass 
ratio binaries in the frequency band of LISA. This may indicate that the 22PN 
waveforms will lead to the accuracy of the data analysis of EMRIs comparable to 
the one resulting from numerical waveforms. 

4.3. Hybrid formula of the energy flux to infinity for a Kerr black hole 

To investigate how higher post-Newtonian order terms for a non-spinning black 
hole improve the accuracy of the 4PN formula for a spinning black hole, 44 ) we con- 
struct a hybrid formula for the energy flux in the post-Newtonian approximation by 
supplementing the 4PN energy flux for a Kerr black hole with our 22PN energy flux 
for a Schwarzschild black hole. In Figs. 3 and 4, ra-PN energy flux includes n-PN 
energy flux to infinity for a Schwarzschild black hole and spin dependent terms of 
the 4PN expression for the energy flux to infinity for a Kerr black hole. The numer- 
ical flux for a Kerr black hole is computed using the same technique 39 )' 40 ) used in 
Sec. 4.1. For the numerical calculation, we set £ = 30 which gives the relative error 
better than 10 -5 up to ISCO for the parameters q used in Figs. 3 and 4. 

In Fig. 3 (Fig. 4), the comparisons for q > {q < 0) are shown. From these 
figures, one will find that for |g| < 0.1 adding higher post-Newtonian order terms for 
a non-spinning black hole achieves the accuracy that is about an order of magnitude 
better than the one using 4PN expression of the energy flux for a spinning black hole. 
For |g| > 0.1, however, the improvement saturates when the PN order is higher than 
8PN. For q < —0.1 the improvement becomes small when \q\ large. For q = —0.9 
one finds at most a factor of two improvement by adding the higher order terms for 
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2 4 6 8 10 12 14 16 18 20 22 24 2 4 6 8 10 12 14 16 18 20 22 24 

t[Month] t[Month] 

Fig. 2. Absolute values of the dephasing during two-year inspiral between the PN and the numerical 
waveforms for the dominant I = m = 2 mode as a function of time in month. The left (right) 
panel represents the early (late) inspiral phase in the LISA band. The left panel shows the 
dephase for System-I of masses (M, fx) = (10 5 , 1Q)Mq, which evolves from ro — 29M to ro — 
16M with associated frequencies /gw £ [4 x 10~ 3 , 10~ 2 ]Hz. The right panel shows the dephase 
for System-II of masses (M, fi) = (10 6 , 10)Mq, which explores orbital radius in a region ro/M £ 
[6.0, 11] and frequencies in a range /gw £ [1-8 x 10 _3 ,4.4 x 10 _3 ]Hz. Note that the dephase 
between the 18PN (22PN) waveforms and numerical waveforms for System-I due to the two-year 
inspiral is about 8 x 10~ 9 (10 -9 ) rads, which is below the lowest value of the dephase in the 
left panel. The dephase between the 22PN waveforms and numerical waveforms for System-II 
after the two-year inspiral is about 10~ 2 rads, which may suggest that the 22PN waveforms will 
provide the data analysis accuracy comparable to the one using numerical waveforms. 



a non-spinning black hole. For q = 0.05, the relative errors of the energy flux around 
ISCO are 2 x 10~\ 4 x 10~ 2 and 6 x 10" 3 for 4PN, 5.5PN and 8PN respectively. 
For q = 0.05, adding 8.5PN or higher order terms does not achieve better accuracy 
than adding 8PN terms, but achieves better accuracy than adding 5.5PN terms. For 
q > 0.1, adding higher post-Newtonian order terms for a non-spinning black hole 
does not always improve the accuracy. For q = 0.1, the relative errors of the 6PN 
energy flux around ISCO becomes larger than that of the 5.5PN energy flux around 
ISCO, but smaller than that of the 4PN energy flux around ISCO. For q = 0.5, the 
relative error of the 6PN energy flux around ISCO is almost same as that of the 4PN 
energy flux around ISCO. For q = 0.9, the relative error of the 6PN (5.5PN) energy 
flux around ISCO is two times larger (an order of magnitude smaller) than that of 
the 4PN energy flux. 

§5. Summary 

Using a systematic method to compute homogeneous solutions of the Teukolsky 
equation, 32 )' 33 ^ we have derived the 22PN expressions of gravitational waves for a 
test particle in circular orbits around a Schwarzschild black hole. The comparison 
of the energy flux between the PN expansion and high precision numerical results 
shows that the relative error of the 22PN energy flux around ISCO is about 10~ 5 . 
We also find that the dephase between the 22PN waveforms and the numerical 
waveforms after two-year inspiral will be smaller than 10~ 2 rads for the most of 
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Fig. 3. Absolute values of the difference of energy flux to infinity between numerical results and PN 
approximation as a function of orbital velocity, v = (M/ro) 1//2 [l + q (M/ro) 3/,2 ] _1//3 , up to ISCO 
for q = 0.01, 0.05, 0.1, 0.3, 0.5 and 0.9. The energy flux at n-PN combines n-PN energy flux for 
a Schwarzschild black hole and spin dependent terms of the 4PN energy flux for a Kerr black 
hole. 44 ' For q < 0.1, higher PN order terms for a non-spinning black hole improve the relative 
accuracy of the energy flux. For q > 0.1, however, higher PN order terms for a non-spinning 
black hole do not always improve the accuracy although one can find the accuracy of some PN 
order is an order of magnitude better than that of the 4PN energy flux. 



the parameter space of EMRIs. This may imply that the 22PN waveforms will 
provide the accuracy in LISA-type data analysis comparable to the one using high 
precision numerical waveforms. (See Ref. 27) for the same discussion using factorized 
resummed waveforms, developed in Ref. 28), at 14PN.) 

The next application will be the case in which a particle moves in circular or- 
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Fig. 4. Same as Fig. 3 but for q = -0.01, -0.05, -0.1, -0.3, -0.5 and q = -0.9. For \q\ < 0.1, 
adding higher post-Newtonian order terms for a non-spinning black hole achieves about an order 
of magnitude better accuracy than using the 4PN energy flux for a spinning black hole. For 
q < —0.1 the improvement becomes small when |<j| large. 



bits around a Kerr black hole. The current available result for the case is 4PN, 44 ) 
which gives the relative errors of 10" 1 (O(l)) for q = 0.1 (q = 0.9) around ISCO 
(See Figs. 3 and 4). To investigate whether our 22PN formula of the energy flux 
improves the accuracy for the case of a Kerr black hole, we construct a hybrid for- 
mula of the energy flux by supplementing the 4PN energy flux for a Kerr black 
hole with the 22PN energy flux for a Schwarzschild black hole. We found that for 
q < 0.1 non-spinning terms of higher post-Newtonian order expressions improve the 
accuracy around ISCO. However, since for q > 0.1 non-spinning terms of higher 
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post-Newtonian order expressions of the total energy flux do not always improve 
the accuracy for a large value of the spin, it is necessary to obtain high PN order 
expressions for the case of a Kerr black hole. For q = 0.9, the numerical calculation 
shows that one needs to include I = 30 modes to achieve the relative error of the 
total energy flux as 10~ 5 around ISCO (See Sec. 4.3). This may indicate that we 
need to compute at least up to 28PN to obtain the relative error of 10 -5 around 
ISCO for q = 0.9. It may be difficult to perform such a high post-Newtonian order 
calculation using our current code since the number of terms necessary to derive 
the PN expressions grows exponentially when the PN order becomes higher (See 
Sec. 1 and Ref. 27)). If one can not obtain sufficiently high post-Newtonian order 
expressions, another approach may be the effective-one-body formalism which can 
include unknown coefficients in the post-Newtonian expressions by calibrating them 
with numerical results. 41 )' 45 ^ 

Acknowledgements 

It is a pleasure to thank Bala R. Iyer for continuous encouragement to look into 
this problem and useful comments on the manuscript. The author is grateful for the 
support of the European Union FEDER funds, the Spanish Ministry of Economy and 
Competitiveness (projects FPA2010-16495 and CSD2007-00042) and the Conselleria 
d'Economia Hisenda i Innovacio of the Govern de les Illes Balears. 

Appendix A 

7PN formula for the energy flux to infinity 



The total energy flux to infinity at 7PN becomes as 
1247 o , 44711 A 8191 



dE fdE\ 



dt V dt 



N 



1 v 2 + 4 ttv 3 v 4 ttv 5 

336 9072 672 



f 6643739519 1712 3424, 16 , 1712, , J , 16285 7 

+ <^ 7 ln(2) H vr 2 ln(v) } v 6 ttv 7 
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177293 , , A n 
1176 V ' ) 

f 256 4 37744140625 2067586193789233570693 

+ \~ 15" 71 " 260941824 n ' ' + 602387400044430000 

246137536815857 27392 437114506833 

157329572400 7 ~ 105 ^ ' 789268480 ^ ' 

271272899815409 , 5861888 , 54784 , . . 2 

ln(2) H ln(2) 7 ln(2) vr 2 

157329572400 v ; 11025 w ' 315 w 

3803225263 2 27392 2 5861888 , . , 2 1465472 2 

H vr 2 7T 7 H ln(2) 2 H 7 2 

10478160 315 ' 11025 v ; 11025 ' 
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where 7 is the Euler constant and C( n ) is the Zeta function. 

We note that coefficients in the post-Newtonian expansion for the total energy 
flux are polynomial functions of tt, 7, ln(n), ln(f) and C( n )- One will find that the 
coefficient at 1.5PN, Att, is derived from the PN expansion of e~ nt / 2 in Eq. (3- 11a) 
for I = m = 2 mode. One will also find that 7, ln(2) and ln(v) terms at 3PN and 
C(3) term at 6PN are derived by the PN expansion of the homogeneous Teukolsky 
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solution in a series of Coulomb wave functions Eq. (3-1) for I = m = 2 mode. 
Similarly, one will find that ln(n) and C(2n ~ 1) terms appear from (n + 1)-PN and 
3n-PN respectively (See Ref. 29)). We also note that (ln(f)) n terms appear from 
3n-PN (See Appendix B and Ref. 29)). One of the reasons is that (ln(v)) n terms 
at 3n-PN are produced by the PN expansion of z v in the homogeneous Teukolsky 
solution in a series of Coulomb wave functions Eq. (3-1), where z = ur = 0(v) and 
v = £ + v^ 1 v G + 0(v 9 ) . Noting that the energy flux is computed by the square of the 
homogeneous Teukolsky solution and = —856/105 for the dominant mode £ = 
m = 2, 31 )~ 33 ) one can estimate leading (ln(u)) n terms in the energy flux using series 
expansion of v ~ 17l2 / 105v an d find that these (ln(f)) n terms at 3n-PN agree with 
the ones in our 22PN energy flux. The explicit expressions for these leading (ln(v)) n 
terms at 3n-PN are already derived in Eq. (44) of Ref. 46) using the renormalization 
group equations and agree with the ones in our 22PN energy flux. 



The total energy flux to infinity at 22PN, which can be used for numerical 
calculation for double precision calculation, becomes as 



1 - 3.7113095238095238 v 2 + 12.566370614359173 v 3 
-4.9284611992945326 v 4 - 38.292835454693446 v 5 
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+ {-7781.9735576745461 +2160.7196071918322 hi(v)}v 1: 
+ {15384.154335011960 + 267.42096713731240 ln(u) 
-340.30387431162941 \n(v) 2 }v 14 
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22PN energy flux to infinity for numerical calculation 
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+52809.58085577387 ln(v) 2 + 113193.57722837660 ln(v) 3 } v 25 
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